home *** CD-ROM | disk | FTP | other *** search
/ io Programmo 60 / IOPROG_60.ISO / soft / c++ / gsl-1.1.1-setup.exe / {app} / src / specfunc / trig.c < prev    next >
Encoding:
C/C++ Source or Header  |  2002-04-18  |  19.4 KB  |  759 lines
  1. /* specfunc/trig.c
  2.  * 
  3.  * Copyright (C) 1996, 1997, 1998, 1999, 2000 Gerard Jungman
  4.  * 
  5.  * This program is free software; you can redistribute it and/or modify
  6.  * it under the terms of the GNU General Public License as published by
  7.  * the Free Software Foundation; either version 2 of the License, or (at
  8.  * your option) any later version.
  9.  * 
  10.  * This program is distributed in the hope that it will be useful, but
  11.  * WITHOUT ANY WARRANTY; without even the implied warranty of
  12.  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  13.  * General Public License for more details.
  14.  * 
  15.  * You should have received a copy of the GNU General Public License
  16.  * along with this program; if not, write to the Free Software
  17.  * Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
  18.  */
  19.  
  20. /* Author:  G. Jungman */
  21.  
  22. #include <config.h>
  23. #include <gsl/gsl_math.h>
  24. #include <gsl/gsl_errno.h>
  25. #include <gsl/gsl_sf_log.h>
  26. #include <gsl/gsl_sf_trig.h>
  27.  
  28. #include "error.h"
  29.  
  30. #include "chebyshev.h"
  31. #include "cheb_eval.c"
  32.  
  33. /* sinh(x) series
  34.  * double-precision for |x| < 1.0
  35.  */
  36. inline
  37. static
  38. int
  39. sinh_series(const double x, double * result)
  40. {
  41.   const double y = x*x;
  42.   const double c0 = 1.0/6.0;
  43.   const double c1 = 1.0/120.0;
  44.   const double c2 = 1.0/5040.0;
  45.   const double c3 = 1.0/362880.0;
  46.   const double c4 = 1.0/39916800.0;
  47.   const double c5 = 1.0/6227020800.0;
  48.   const double c6 = 1.0/1307674368000.0;
  49.   const double c7 = 1.0/355687428096000.0;
  50.   *result = x*(1.0 + y*(c0+y*(c1+y*(c2+y*(c3+y*(c4+y*(c5+y*(c6+y*c7))))))));
  51.   return GSL_SUCCESS;
  52. }
  53.  
  54.  
  55. /* cosh(x)-1 series
  56.  * double-precision for |x| < 1.0
  57.  */
  58. inline
  59. static
  60. int
  61. cosh_m1_series(const double x, double * result)
  62. {
  63.   const double y = x*x;
  64.   const double c0 = 0.5;
  65.   const double c1 = 1.0/24.0;
  66.   const double c2 = 1.0/720.0;
  67.   const double c3 = 1.0/40320.0;
  68.   const double c4 = 1.0/3628800.0;
  69.   const double c5 = 1.0/479001600.0;
  70.   const double c6 = 1.0/87178291200.0;
  71.   const double c7 = 1.0/20922789888000.0;
  72.   const double c8 = 1.0/6402373705728000.0;
  73.   *result = y*(c0+y*(c1+y*(c2+y*(c3+y*(c4+y*(c5+y*(c6+y*(c7+y*c8))))))));
  74.   return GSL_SUCCESS;
  75. }
  76.  
  77.  
  78. /* Chebyshev expansion for f(t) = sinc((t+1)/2), -1 < t < 1
  79.  */
  80. static double sinc_data[17] = {
  81.   1.133648177811747875422,
  82.  -0.532677564732557348781,
  83.  -0.068293048346633177859,
  84.   0.033403684226353715020,
  85.   0.001485679893925747818,
  86.  -0.000734421305768455295,
  87.  -0.000016837282388837229,
  88.   0.000008359950146618018,
  89.   0.000000117382095601192,
  90.  -0.000000058413665922724,
  91.  -0.000000000554763755743,
  92.   0.000000000276434190426,
  93.   0.000000000001895374892,
  94.  -0.000000000000945237101,
  95.  -0.000000000000004900690,
  96.   0.000000000000002445383,
  97.   0.000000000000000009925
  98. };
  99. static cheb_series sinc_cs = {
  100.   sinc_data,
  101.   16,
  102.   -1, 1,
  103.   10
  104. };
  105.  
  106.  
  107. /* Chebyshev expansion for f(t) = g((t+1)Pi/8), -1<t<1
  108.  * g(x) = (sin(x)/x - 1)/(x*x)
  109.  */
  110. static double sin_data[12] = {
  111.   -0.3295190160663511504173,
  112.    0.0025374284671667991990,
  113.    0.0006261928782647355874,
  114.   -4.6495547521854042157541e-06,
  115.   -5.6917531549379706526677e-07,
  116.    3.7283335140973803627866e-09,
  117.    3.0267376484747473727186e-10,
  118.   -1.7400875016436622322022e-12,
  119.   -1.0554678305790849834462e-13,
  120.    5.3701981409132410797062e-16,
  121.    2.5984137983099020336115e-17,
  122.   -1.1821555255364833468288e-19
  123. };
  124. static cheb_series sin_cs = {
  125.   sin_data,
  126.   11,
  127.   -1, 1,
  128.   11
  129. };
  130.  
  131. /* Chebyshev expansion for f(t) = g((t+1)Pi/8), -1<t<1
  132.  * g(x) = (2(cos(x) - 1)/(x^2) + 1) / x^2
  133.  */
  134. static double cos_data[11] = {
  135.   0.165391825637921473505668118136,
  136.  -0.00084852883845000173671196530195,
  137.  -0.000210086507222940730213625768083,
  138.   1.16582269619760204299639757584e-6,
  139.   1.43319375856259870334412701165e-7,
  140.  -7.4770883429007141617951330184e-10,
  141.  -6.0969994944584252706997438007e-11,
  142.   2.90748249201909353949854872638e-13,
  143.   1.77126739876261435667156490461e-14,
  144.  -7.6896421502815579078577263149e-17,
  145.  -3.7363121133079412079201377318e-18
  146. };
  147. static cheb_series cos_cs = {
  148.   cos_data,
  149.   10,
  150.   -1, 1,
  151.   10
  152. };
  153.  
  154.  
  155. /*-*-*-*-*-*-*-*-*-*-*-* Functions with Error Codes *-*-*-*-*-*-*-*-*-*-*-*/
  156.  
  157. /* I would have prefered just using the library sin() function.
  158.  * But after some experimentation I decided that there was
  159.  * no good way to understand the error; library sin() is just a black box.
  160.  * So we have to roll our own.
  161.  */
  162. int
  163. gsl_sf_sin_e(double x, gsl_sf_result * result)
  164. {
  165.   /* CHECK_POINTER(result) */
  166.  
  167.   {
  168.     const double P1 = 7.85398125648498535156e-1;
  169.     const double P2 = 3.77489470793079817668e-8;
  170.     const double P3 = 2.69515142907905952645e-15;
  171.  
  172.     const double sgn_x = GSL_SIGN(x);
  173.     const double abs_x = fabs(x);
  174.  
  175.     if(abs_x < GSL_ROOT4_DBL_EPSILON) {
  176.       const double x2 = x*x;
  177.       result->val = x * (1.0 - x2/6.0);
  178.       result->err = fabs(x*x2*x2 / 100.0);
  179.       return GSL_SUCCESS;
  180.     }
  181.     else {
  182.       double sgn_result = sgn_x;
  183.       double y = floor(abs_x/(0.25*M_PI));
  184.       int octant = y - ldexp(floor(ldexp(y,-3)),3);
  185.       int stat_cs;
  186.       double z;
  187.  
  188.       if(GSL_IS_ODD(octant)) {
  189.         octant += 1;
  190.     octant &= 07;
  191.     y += 1.0;
  192.       }
  193.  
  194.       if(octant > 3) {
  195.         octant -= 4;
  196.     sgn_result = -sgn_result;
  197.       }
  198.       
  199.       z = ((abs_x - y * P1) - y * P2) - y * P3;
  200.  
  201.       if(octant == 0) {
  202.         gsl_sf_result sin_cs_result;
  203.     const double t = 8.0*fabs(z)/M_PI - 1.0;
  204.     stat_cs = cheb_eval_e(&sin_cs, t, &sin_cs_result);
  205.         result->val = z * (1.0 + z*z * sin_cs_result.val);
  206.       }
  207.       else { /* octant == 2 */
  208.         gsl_sf_result cos_cs_result;
  209.     const double t = 8.0*fabs(z)/M_PI - 1.0;
  210.     stat_cs = cheb_eval_e(&cos_cs, t, &cos_cs_result);
  211.         result->val = 1.0 - 0.5*z*z * (1.0 - z*z * cos_cs_result.val);
  212.       }
  213.  
  214.       result->val *= sgn_result;
  215.  
  216.       if(abs_x > 1.0/GSL_DBL_EPSILON) {
  217.         result->err = fabs(result->val);
  218.       }
  219.       else if(abs_x > 100.0/GSL_SQRT_DBL_EPSILON) {
  220.         result->err = 2.0 * abs_x * GSL_DBL_EPSILON * fabs(result->val);
  221.       }
  222.       else if(abs_x > 0.1/GSL_SQRT_DBL_EPSILON) {
  223.         result->err = 2.0 * GSL_SQRT_DBL_EPSILON * fabs(result->val);
  224.       }
  225.       else {
  226.         result->err = 2.0 * GSL_DBL_EPSILON * fabs(result->val);
  227.       }
  228.  
  229.       return stat_cs;
  230.     }
  231.   }
  232. }
  233.  
  234.  
  235. int
  236. gsl_sf_cos_e(double x, gsl_sf_result * result)
  237. {
  238.   /* CHECK_POINTER(result) */
  239.  
  240.   {
  241.     const double P1 = 7.85398125648498535156e-1;
  242.     const double P2 = 3.77489470793079817668e-8;
  243.     const double P3 = 2.69515142907905952645e-15;
  244.  
  245.     const double abs_x = fabs(x);
  246.  
  247.     if(abs_x < GSL_ROOT4_DBL_EPSILON) {
  248.       const double x2 = x*x;
  249.       result->val = 1.0 - 0.5*x2;
  250.       result->err = fabs(x2*x2/12.0);
  251.       return GSL_SUCCESS;
  252.     }
  253.     else {
  254.       double sgn_result = 1.0;
  255.       double y = floor(abs_x/(0.25*M_PI));
  256.       int octant = y - ldexp(floor(ldexp(y,-3)),3);
  257.       int stat_cs;
  258.       double z;
  259.  
  260.       if(GSL_IS_ODD(octant)) {
  261.         octant += 1;
  262.     octant &= 07;
  263.     y += 1.0;
  264.       }
  265.  
  266.       if(octant > 3) {
  267.         octant -= 4;
  268.     sgn_result = -sgn_result;
  269.       }
  270.  
  271.       if(octant > 1) {
  272.         sgn_result = -sgn_result;
  273.       }
  274.  
  275.       z = ((abs_x - y * P1) - y * P2) - y * P3;
  276.  
  277.       if(octant == 0) {
  278.         gsl_sf_result cos_cs_result;
  279.         const double t = 8.0*fabs(z)/M_PI - 1.0;
  280.         stat_cs = cheb_eval_e(&cos_cs, t, &cos_cs_result);
  281.         result->val = 1.0 - 0.5*z*z * (1.0 - z*z * cos_cs_result.val);
  282.       }
  283.       else { /* octant == 2 */
  284.         gsl_sf_result sin_cs_result;
  285.         const double t = 8.0*fabs(z)/M_PI - 1.0;
  286.         stat_cs = cheb_eval_e(&sin_cs, t, &sin_cs_result);
  287.         result->val = z * (1.0 + z*z * sin_cs_result.val);
  288.       }
  289.  
  290.       result->val *= sgn_result;
  291.  
  292.       if(abs_x > 1.0/GSL_DBL_EPSILON) {
  293.         result->err = fabs(result->val);
  294.       }
  295.       else if(abs_x > 100.0/GSL_SQRT_DBL_EPSILON) {
  296.         result->err = 2.0 * abs_x * GSL_DBL_EPSILON * fabs(result->val);
  297.       }
  298.       else if(abs_x > 0.1/GSL_SQRT_DBL_EPSILON) {
  299.         result->err = 2.0 * GSL_SQRT_DBL_EPSILON * fabs(result->val);
  300.       }
  301.       else {
  302.         result->err = 2.0 * GSL_DBL_EPSILON * fabs(result->val);
  303.       }
  304.  
  305.       return stat_cs;
  306.     }
  307.   }
  308. }
  309.  
  310.  
  311. int
  312. gsl_sf_hypot_e(const double x, const double y, gsl_sf_result * result)
  313. {
  314.   /* CHECK_POINTER(result) */
  315.  
  316.   if(x == 0.0 && y == 0.0) {
  317.     result->val = 0.0;
  318.     result->err = 0.0;
  319.     return GSL_SUCCESS;
  320.   }
  321.   else {
  322.     const double a = fabs(x);
  323.     const double b = fabs(y);
  324.     const double min = GSL_MIN_DBL(a,b);
  325.     const double max = GSL_MAX_DBL(a,b);
  326.     const double rat = min/max;
  327.     const double root_term = sqrt(1.0 + rat*rat);
  328.  
  329.     if(max < GSL_DBL_MAX/root_term) {
  330.       result->val = max * root_term;
  331.       result->err = 2.0 * GSL_DBL_EPSILON * fabs(result->val);
  332.       return GSL_SUCCESS;
  333.     }
  334.     else {
  335.       OVERFLOW_ERROR(result);
  336.     }
  337.   }
  338. }
  339.  
  340.  
  341. int
  342. gsl_sf_complex_sin_e(const double zr, const double zi,
  343.                         gsl_sf_result * szr, gsl_sf_result * szi)
  344. {
  345.   /* CHECK_POINTER(szr) */
  346.   /* CHECK_POINTER(szi) */
  347.  
  348.   if(fabs(zi) < 1.0) {
  349.     double ch_m1, sh;
  350.     sinh_series(zi, &sh);
  351.     cosh_m1_series(zi, &ch_m1);
  352.     szr->val = sin(zr)*(ch_m1 + 1.0);
  353.     szi->val = cos(zr)*sh;
  354.     szr->err = 2.0 * GSL_DBL_EPSILON * fabs(szr->val);
  355.     szi->err = 2.0 * GSL_DBL_EPSILON * fabs(szi->val);
  356.     return GSL_SUCCESS;
  357.   }
  358.   else if(fabs(zi) < GSL_LOG_DBL_MAX) {
  359.     double ex = exp(zi);
  360.     double ch = 0.5*(ex+1.0/ex);
  361.     double sh = 0.5*(ex-1.0/ex);
  362.     szr->val = sin(zr)*ch;
  363.     szi->val = cos(zr)*sh;
  364.     szr->err = 2.0 * GSL_DBL_EPSILON * fabs(szr->val);
  365.     szi->err = 2.0 * GSL_DBL_EPSILON * fabs(szi->val);
  366.     return GSL_SUCCESS;
  367.   }
  368.   else {
  369.     OVERFLOW_ERROR_2(szr, szi);
  370.   }
  371. }
  372.  
  373.  
  374. int
  375. gsl_sf_complex_cos_e(const double zr, const double zi,
  376.                         gsl_sf_result * czr, gsl_sf_result * czi)
  377. {
  378.   /* CHECK_POINTER(czr) */
  379.   /* CHECK_POINTER(czi) */
  380.  
  381.   if(fabs(zi) < 1.0) {
  382.     double ch_m1, sh;
  383.     sinh_series(zi, &sh);
  384.     cosh_m1_series(zi, &ch_m1);
  385.     czr->val =  cos(zr)*(ch_m1 + 1.0);
  386.     czi->val = -sin(zr)*sh;
  387.     czr->err = 2.0 * GSL_DBL_EPSILON * fabs(czr->val);
  388.     czi->err = 2.0 * GSL_DBL_EPSILON * fabs(czi->val);
  389.     return GSL_SUCCESS;
  390.   }
  391.   else if(fabs(zi) < GSL_LOG_DBL_MAX) {
  392.     double ex = exp(zi);
  393.     double ch = 0.5*(ex+1.0/ex);
  394.     double sh = 0.5*(ex-1.0/ex);
  395.     czr->val =  cos(zr)*ch;
  396.     czi->val = -sin(zr)*sh;
  397.     czr->err = 2.0 * GSL_DBL_EPSILON * fabs(czr->val);
  398.     czi->err = 2.0 * GSL_DBL_EPSILON * fabs(czi->val);
  399.     return GSL_SUCCESS;
  400.   }
  401.   else {
  402.     OVERFLOW_ERROR_2(czr,czi);
  403.   }
  404. }
  405.  
  406.  
  407. int
  408. gsl_sf_complex_logsin_e(const double zr, const double zi,
  409.                            gsl_sf_result * lszr, gsl_sf_result * lszi)
  410. {
  411.   /* CHECK_POINTER(lszr) */
  412.   /* CHECK_POINTER(lszi) */
  413.  
  414.   if(zi > 60.0) {
  415.     lszr->val = -M_LN2 + zi;
  416.     lszi->val =  0.5*M_PI - zr;
  417.     lszr->err = 2.0 * GSL_DBL_EPSILON * fabs(lszr->val);
  418.     lszi->err = 2.0 * GSL_DBL_EPSILON * fabs(lszi->val);
  419.   }
  420.   else if(zi < -60.0) {
  421.     lszr->val = -M_LN2 - zi;
  422.     lszi->val = -0.5*M_PI + zr;
  423.     lszr->err = 2.0 * GSL_DBL_EPSILON * fabs(lszr->val);
  424.     lszi->err = 2.0 * GSL_DBL_EPSILON * fabs(lszi->val);
  425.   }
  426.   else {
  427.     gsl_sf_result sin_r, sin_i;
  428.     int status;
  429.     gsl_sf_complex_sin_e(zr, zi, &sin_r, &sin_i); /* ok by construction */
  430.     status = gsl_sf_complex_log_e(sin_r.val, sin_i.val, lszr, lszi);
  431.     if(status == GSL_EDOM) {
  432.       DOMAIN_ERROR_2(lszr, lszi);
  433.     }
  434.   }
  435.   return gsl_sf_angle_restrict_symm_e(&(lszi->val));
  436. }
  437.  
  438.  
  439. int
  440. gsl_sf_lnsinh_e(const double x, gsl_sf_result * result)
  441. {
  442.   /* CHECK_POINTER(result) */
  443.  
  444.   if(x <= 0.0) {
  445.     DOMAIN_ERROR(result);
  446.   }
  447.   else if(fabs(x) < 1.0) {
  448.     double eps;
  449.     sinh_series(x, &eps);
  450.     result->val = log(eps);
  451.     result->err = 2.0 * GSL_DBL_EPSILON * fabs(result->val);
  452.     return GSL_SUCCESS;
  453.   }
  454.   else if(x < -0.5*GSL_LOG_DBL_EPSILON) {
  455.     result->val = x + log(0.5*(1.0 - exp(-2.0*x)));
  456.     result->err = 2.0 * GSL_DBL_EPSILON * fabs(result->val);
  457.     return GSL_SUCCESS;
  458.   }
  459.   else {
  460.     result->val = -M_LN2 + x;
  461.     result->err = 2.0 * GSL_DBL_EPSILON * fabs(result->val);
  462.     return GSL_SUCCESS;
  463.   }
  464. }
  465.  
  466.  
  467. int gsl_sf_lncosh_e(const double x, gsl_sf_result * result)
  468. {
  469.   /* CHECK_POINTER(result) */
  470.  
  471.   if(fabs(x) < 1.0) {
  472.     double eps;
  473.     cosh_m1_series(x, &eps);
  474.     return gsl_sf_log_1plusx_e(eps, result);
  475.   }
  476.   else if(x < -0.5*GSL_LOG_DBL_EPSILON) {
  477.     result->val = x + log(0.5*(1.0 + exp(-2.0*x)));
  478.     result->err = 2.0 * GSL_DBL_EPSILON * fabs(result->val);
  479.     return GSL_SUCCESS;
  480.   }
  481.   else {
  482.     result->val = -M_LN2 + x;
  483.     result->err = 2.0 * GSL_DBL_EPSILON * fabs(result->val);
  484.     return GSL_SUCCESS;
  485.   }
  486. }
  487.  
  488.  
  489. /*
  490. inline int gsl_sf_sincos_e(const double theta, double * s, double * c)
  491. {
  492.   double tan_half = tan(0.5 * theta);
  493.   double den = 1. + tan_half*tan_half;
  494.   double cos_theta = (1.0 - tan_half*tan_half) / den;
  495.   double sin_theta = 2.0 * tan_half / den;
  496. }
  497. */
  498.  
  499. int
  500. gsl_sf_polar_to_rect(const double r, const double theta,
  501.                           gsl_sf_result * x, gsl_sf_result * y)
  502. {
  503.   double t   = theta;
  504.   int status = gsl_sf_angle_restrict_symm_e(&t);
  505.   double c = cos(t);
  506.   double s = sin(t);
  507.   x->val = r * cos(t);
  508.   y->val = r * sin(t);
  509.   x->err  = r * fabs(s * GSL_DBL_EPSILON * t);
  510.   x->err += 2.0 * GSL_DBL_EPSILON * fabs(x->val);
  511.   y->err  = r * fabs(c * GSL_DBL_EPSILON * t);
  512.   y->err += 2.0 * GSL_DBL_EPSILON * fabs(y->val);
  513.   return status;
  514. }
  515.  
  516.  
  517. int
  518. gsl_sf_rect_to_polar(const double x, const double y,
  519.                           gsl_sf_result * r, gsl_sf_result * theta)
  520. {
  521.   int stat_h = gsl_sf_hypot_e(x, y, r);
  522.   if(r->val > 0.0) {
  523.     theta->val = atan2(y, x);
  524.     theta->err = 2.0 * GSL_DBL_EPSILON * fabs(theta->val);
  525.     return stat_h;
  526.   }
  527.   else {
  528.     DOMAIN_ERROR(theta);
  529.   }
  530. }
  531.  
  532.  
  533. int gsl_sf_angle_restrict_symm_err_e(const double theta, gsl_sf_result * result)
  534. {
  535.   /* synthetic extended precision constants */
  536.   const double P1 = 4 * 7.8539812564849853515625e-01;
  537.   const double P2 = 4 * 3.7748947079307981766760e-08;
  538.   const double P3 = 4 * 2.6951514290790594840552e-15;
  539.   const double TwoPi = 2*(P1 + P2 + P3);
  540.  
  541.   const double y = 2*floor(theta/TwoPi);
  542.   double r = ((theta - y*P1) - y*P2) - y*P3;
  543.  
  544.   if(r >  M_PI) r -= TwoPi;
  545.   result->val = r;
  546.  
  547.   if(theta > 0.0625/GSL_DBL_EPSILON) {
  548.     result->err = fabs(result->val);
  549.     GSL_ERROR ("error", GSL_ELOSS);
  550.   }
  551.   else if(theta > 0.0625/GSL_SQRT_DBL_EPSILON) {
  552.     result->err = GSL_SQRT_DBL_EPSILON * fabs(result->val);
  553.     return GSL_SUCCESS;
  554.   }
  555.   else {
  556.     result->err = 2.0 * GSL_DBL_EPSILON * fabs(result->val);
  557.     return GSL_SUCCESS;
  558.   }
  559. }
  560.  
  561.  
  562. int gsl_sf_angle_restrict_pos_err_e(const double theta, gsl_sf_result * result)
  563. {
  564.   /* synthetic extended precision constants */
  565.   const double P1 = 4 * 7.85398125648498535156e-01;
  566.   const double P2 = 4 * 3.77489470793079817668e-08;
  567.   const double P3 = 4 * 2.69515142907905952645e-15;
  568.   const double TwoPi = 2*(P1 + P2 + P3);
  569.  
  570.   const double y = 2*floor(theta/TwoPi);
  571.  
  572.   result->val = ((theta - y*P1) - y*P2) - y*P3;
  573.  
  574.   if(theta > 0.0625/GSL_DBL_EPSILON) {
  575.     result->err = fabs(result->val);
  576.     GSL_ERROR ("error", GSL_ELOSS);
  577.   }
  578.   else if(theta > 0.0625/GSL_SQRT_DBL_EPSILON) {
  579.     result->err = GSL_SQRT_DBL_EPSILON * fabs(result->val);
  580.     return GSL_SUCCESS;
  581.   }
  582.   else {
  583.     result->err = 2.0 * GSL_DBL_EPSILON * fabs(result->val);
  584.     return GSL_SUCCESS;
  585.   }
  586. }
  587.  
  588.  
  589. int gsl_sf_angle_restrict_symm_e(double * theta)
  590. {
  591.   gsl_sf_result r;
  592.   int stat = gsl_sf_angle_restrict_symm_err_e(*theta, &r);
  593.   *theta = r.val;
  594.   return stat;
  595. }
  596.  
  597.  
  598. int gsl_sf_angle_restrict_pos_e(double * theta)
  599. {
  600.   gsl_sf_result r;
  601.   int stat = gsl_sf_angle_restrict_pos_err_e(*theta, &r);
  602.   *theta = r.val;
  603.   return stat;
  604. }
  605.  
  606.  
  607. int gsl_sf_sin_err_e(const double x, const double dx, gsl_sf_result * result)
  608. {
  609.   int stat_s = gsl_sf_sin_e(x, result);
  610.   result->err += fabs(cos(x) * dx);
  611.   result->err += GSL_DBL_EPSILON * fabs(result->val);
  612.   return stat_s;
  613. }
  614.  
  615.  
  616. int gsl_sf_cos_err_e(const double x, const double dx, gsl_sf_result * result)
  617. {
  618.   int stat_c = gsl_sf_cos_e(x, result);
  619.   result->err += fabs(sin(x) * dx);
  620.   result->err += GSL_DBL_EPSILON * fabs(result->val);
  621.   return stat_c;
  622. }
  623.  
  624.  
  625. #if 0
  626. int
  627. gsl_sf_sin_pi_x_e(const double x, gsl_sf_result * result)
  628. {
  629.   /* CHECK_POINTER(result) */
  630.  
  631.   if(-100.0 < x && x < 100.0) {
  632.     result->val = sin(M_PI * x) / (M_PI * x);
  633.     result->err = 2.0 * GSL_DBL_EPSILON * fabs(result->val);
  634.     return GSL_SUCCESS;
  635.   }
  636.   else {
  637.     const double N = floor(x + 0.5);
  638.     const double f = x - N;
  639.  
  640.     if(N < INT_MAX && N > INT_MIN) {
  641.       /* Make it an integer if we can. Saves another
  642.        * call to floor().
  643.        */
  644.       const int intN    = (int)N;
  645.       const double sign = ( GSL_IS_ODD(intN) ? -1.0 : 1.0 );
  646.       result->val = sign * sin(M_PI * f);
  647.       result->err = GSL_DBL_EPSILON * fabs(result->val);
  648.     }
  649.     else if(N > 2.0/GSL_DBL_EPSILON || N < -2.0/GSL_DBL_EPSILON) {
  650.       /* All integer-valued floating point numbers
  651.        * bigger than 2/eps=2^53 are actually even.
  652.        */
  653.       result->val = 0.0;
  654.       result->err = 0.0;
  655.     }
  656.     else {
  657.       const double resN = N - 2.0*floor(0.5*N); /* 0 for even N, 1 for odd N */
  658.       const double sign = ( fabs(resN) > 0.5 ? -1.0 : 1.0 );
  659.       result->val = sign * sin(M_PI*f);
  660.       result->err = GSL_DBL_EPSILON * fabs(result->val);
  661.     }
  662.  
  663.     return GSL_SUCCESS;
  664.   }
  665. }
  666. #endif
  667.  
  668.  
  669. int gsl_sf_sinc_e(double x, gsl_sf_result * result)
  670. {
  671.   /* CHECK_POINTER(result) */
  672.  
  673.   {
  674.     const double ax = fabs(x);
  675.  
  676.     if(ax < 0.8) {
  677.       /* Do not go to the limit of the fit since
  678.        * there is a zero there and the Chebyshev
  679.        * accuracy will go to zero.
  680.        */
  681.       return cheb_eval_e(&sinc_cs, 2.0*ax-1.0, result);
  682.     }
  683.     else if(ax < 100.0) {
  684.       /* Small arguments are no problem.
  685.        * We trust the library sin() to
  686.        * roughly machine precision.
  687.        */
  688.       result->val = sin(M_PI * ax)/(M_PI * ax);
  689.       result->err = 2.0 * GSL_DBL_EPSILON * fabs(result->val);
  690.       return GSL_SUCCESS;
  691.     }
  692.     else {
  693.       /* Large arguments must be handled separately.
  694.        */
  695.       const double r = M_PI*ax;
  696.       gsl_sf_result s;
  697.       int stat_s = gsl_sf_sin_e(r, &s);
  698.       result->val = s.val/r;
  699.       result->err = s.err/r + 2.0 * GSL_DBL_EPSILON * fabs(result->val);
  700.       return stat_s;
  701.     }
  702.   }
  703. }
  704.  
  705.  
  706.  
  707. /*-*-*-*-*-*-*-*-*-* Functions w/ Natural Prototypes *-*-*-*-*-*-*-*-*-*-*/
  708.  
  709. #include "eval.h"
  710.  
  711. double gsl_sf_sin(const double x)
  712. {
  713.   EVAL_RESULT(gsl_sf_sin_e(x, &result));
  714. }
  715.  
  716. double gsl_sf_cos(const double x)
  717. {
  718.   EVAL_RESULT(gsl_sf_cos_e(x, &result));
  719. }
  720.  
  721. double gsl_sf_hypot(const double x, const double y)
  722. {
  723.   EVAL_RESULT(gsl_sf_hypot_e(x, y, &result));
  724. }
  725.  
  726. double gsl_sf_lnsinh(const double x)
  727. {
  728.   EVAL_RESULT(gsl_sf_lnsinh_e(x, &result));
  729. }
  730.  
  731. double gsl_sf_lncosh(const double x)
  732. {
  733.   EVAL_RESULT(gsl_sf_lncosh_e(x, &result));
  734. }
  735.  
  736. double gsl_sf_angle_restrict_symm(const double theta)
  737. {
  738.   double result = theta;
  739.   EVAL_DOUBLE(gsl_sf_angle_restrict_symm_e(&result));
  740. }
  741.  
  742. double gsl_sf_angle_restrict_pos(const double theta)
  743. {
  744.   double result = theta;
  745.   EVAL_DOUBLE(gsl_sf_angle_restrict_pos_e(&result));
  746. }
  747.  
  748. #if 0
  749. double gsl_sf_sin_pi_x(const double x)
  750. {
  751.   EVAL_RESULT(gsl_sf_sin_pi_x_e(x, &result));
  752. }
  753. #endif
  754.  
  755. double gsl_sf_sinc(const double x)
  756. {
  757.   EVAL_RESULT(gsl_sf_sinc_e(x, &result));
  758. }
  759.